library(sf)
library(dplyr)
library(ggplot2)
library(haven)
library(ggmap)

######### Figure 3 ######### 

# reading data

shape_cty <- st_read("cb_2015_us_county_5m.shp")

blood_center_geo <- read.csv("blood_center_geo.csv")
blood_center_geo_sf <- st_as_sf(blood_center_geo, coords = c("location_lng", "location_lat"), crs = 4326)
blood_center_geo_sf <- st_transform(blood_center_geo_sf, st_crs(shape_cty))

# spatial merging

intersection_state <- st_intersects(shape_cty, blood_center_geo_sf, sparse = FALSE)
blood_center_geo_sf$state_id <- apply(intersection_state, 2, function(x) shape_cty[x,]$STATEFP)

intersection_county <- st_intersects(shape_cty, blood_center_geo_sf, sparse = FALSE)
blood_center_geo_sf$cty_id <- apply(intersection_county, 2, function(x) shape_cty[x,]$COUNTYFP)

# subset blood donation centers

coords <- st_coordinates(blood_center_geo_sf)
longitude <- coords[, 'X']
latitude <- coords[, 'Y']
blood_center_geo_sf$longitude <- longitude
blood_center_geo_sf$latitude <- latitude

map_blood_center_geo2_data <- blood_center_geo_sf %>%
  filter(longitude < -50, longitude > -135, latitude < 50, latitude > 22)

# subset county boundaries

bbox <- st_as_sfc(st_bbox(c(xmin = -135, xmax = -50, ymin = 22, ymax = 50), crs = st_crs(shape_cty)))
shape_cty <- shape_cty[st_intersects(shape_cty, bbox, sparse = FALSE), ]


# map

pdf("Figure3.pdf", width=7.5, height=5)
mapPlot <- ggplot() +
  geom_sf(data = shape_cty, fill = "white", color = "lavenderblush4") +
  geom_sf(data = map_blood_center_geo2_data, aes(geometry = geometry), size = 1 , alpha = 0.5) +
  theme_bw() + 
  theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank() ,axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), legend.position="bottom",    legend.key.width = unit(1, "cm"), legend.text = element_text(size = 8)) +
  coord_sf(crs = "+proj=aea +lat_1=30 +lat_2=60 +lat_0=45 +lon_0=-96")
mapPlot
dev.off()

######### Figure A1 ######### 

anes_aggregate <- read_dta("anes_aggregate.dta")

pdf("FigureA1.pdf",width=5,height=5)

ggplot(data = anes_aggregate, mapping = aes(st_ideo7_var, affective_polarization_avg)) + 
  geom_point() + ylim(0,100) +
  geom_smooth(method="lm") + xlab("Variance of Ideology") + ylab("Affective Polarization")  

dev.off()

